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Abstract 

We present a practical numerical method for evaluating the Lagrange multi- 
pliers necessary for maintaining a constrained linear geometry of particles in 
dynamical simulations. The method involves no iterations, and is limited in 
accuracy only by the numerical methods for solving small systems of linear 
equations. As a result of the non- iterative and exact (within numerical accu- 
racy) nature of the procedure there is no drift in the constrained geometry, and 
the method is therefore readily applied to molecular dynamics simulations of, 
e.g., rigid linear molecules or materials of non-spherical grains. We illustrate 
the approach through implementation in the commonly used second-order 
velocity explicit Verlet method. 
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I. INTRODUCTION 



Computational studies of the dynamics and statistics of large particle ensembles have 
been successfully conducted for several decades [1,2], and many different types of inter- 
particle force fields have been developed for specific physical applications. The particles 
of interest in Molecular Dynamics (MD) or Discrete Element Methods (DEMs) are often 
represented by spatial point coordinates and the associated momenta. For example, atomic- 
scale modeling of materials considers each atom as a point particle, and the nature of the 
material is defined by the interatomic potential. Likewise, the dynamics of inhomogeneous 
atomic ensembles, such as molecular chains in solution, can be modeled by point particles 
with different interaction potentials defining the natiirc of the chemistry between two or 
more atoms. Given a complex ensemble of interactions, the involved time scales of the dy- 
namics may span many orders of magnitude, thereby limiting the efficiency of a simulation 
approach since a useful numerical time step for a simulation is inversely proportional to 
the highest frequency in the system. Thus, if particles interact such that high-frequency, 
small-amplitude oscillations result, it may be advantageous to disregard the dynamics of 
the interaction and constrain the relative degree of freedom between the two objects to 
the mean distance. In molecular systems at room temperature conditions, constraints are 
further justified by the fact that many relevant chemical bond vibrations are improbable 
to be found in an excited quantum mechanical energy state. Other examples of the appli- 
cation of holonomic constraints can be found in studies of granular materials composed of 
non-spherical particles. Many such studies relate to particles of complex geometry or to 
particles that are connected through actual physical constraints. Modeling these systems by 
spherical particles interacting through vibrational potentials is obviously not desirable as it 
may introduce spurious deformations and internal oscillations. 

As a result of the interest in constraints, much work has been devoted in the literature 
to develop practical, accurate, and efficient methods for introducing constraints into the 
numerical methods to study the dynamics of non-spherical particles. Much of the pioneering 
work culminated in the SHAKE [3] and RATTLE [4] algorithms. These methods consider the 
constraints in a (point) particle ensemble as a collection of linearly independent constraints 
on distances between pairs of particles. For any given time step of a numerical temporal 
integrator of the equations of motion for each particle, the constraints are enforced for the 
resulting positions through the self-consistent solution of a set of nonlinear equations (one 
for each constraint). This solution is usually obtained through an iterative approach where 
computational efficiency has to be balanced with the desire for accuracy. These approaches 
have been very successful and are still the core in many applications. Some modifications 
have been introduced for specific systems and geometries. For example, SETTLE [5] provides 
a direct analytical calculation of the three constraints necessary for rigid triangular objects 
(developed for three-point classical water models), and "Fast SHAKE" [6] optimizes the 
iterative method for finding the constraints necessary for modeling small rigid molecules. 
Other suggestions include the EEM approach [7,8], which evaluates the constraints non- 
iteratively at the beginning of a time step where all positions arc known. As a consequence 
of iterative methods and algorithmic approximations, the updated positions do not exactly 
fulfill the holonomic constraints, and corrective measures must be applied if and when the 
cumulative discrepancy becomes unacceptable. 
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Most methods for evaluating constraints have been developed for ensembles of linearly 
independent constraints. Several publications, including [3], consider the linear geometry 
as a separate case due to occurring singularities in the methods for linearly independent 
constraints. A few direct methods for hnear geometry have been outlined (see, e.g., [9]), but 
they are possible only because of perfect three-body symmetry. We demonstrate a simple 
method for evaluating the necessary Lagrange multipliers for constraining any number of 
different particles onto a linear geometry. We show that the resulting method is efficient, 
non-iterative, exact, and self-correcting. 



II. EQUATIONS OF MOTION WITH CONSTRAINTS 

We consider the dynamics of N particles constrained on a straight line (see Figure 1). 
The equations of motion for the coordinate, Rj, of the jth particle is given by: 
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where mj is the mass, friction is represented by aj, and fj is an extra-molecular force 
on the jth particle. Throughout the paper, overbar and overdot "'" denote a vector 
and a temporal derivative, respectively. The two particles, j — 1 and j — N, define the 
direction and length of the object through the scalar Lagrange multiplier Aiat and the vector 
RiN = Rn — Ri- The length Iji is defined by Iji = \Ri — Rj\. All other particles, 1 < j < A^, 
are constrained to their relative positions through the vector multipliers Xji foi i = 1,N, that 
connect each of those particles exclusively to the defining particles of j = 1, A^. Notice that 
all the forces of the constraints cancel for the ensemble. At any given time the constraints 
are enforced by the geometric conditions: 
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where all particles 1 < j < A^ lie between particles 1 and A^, such that Iji + Z^tv = ^lAf- The 
number of geometrical conditions given in (2) and (3) is 3(A^ — 2) -|- 1 (three-dimensional 
space), whereas the number of Lagrange multipliers to be determined for the solution of (1) 
is 6(A^ — 2) -I- 1. Since the constraints should not contribute net torque to the object, we 
further have the conditions for 1 < j < A^: 



RlN X (ijl^jl — IjN^jN^ — 



(4) 



By writing 
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(5) 
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with • RiN — and A'^^ = A'jji?iAr/|-RiAr|, we can express (4) in the form 



(6) 



i.e., the constraint forces, responsible for maintaining particle j at its desired relative posi- 
tion, are balanced in their components orthogonal to the direction of the object such that 
no net torque is contributed. 

However, it is important to realize that this balance does not necessarily apply to the 
longitudinal components, A^j^ and A^^'^y. Since the object is rigid, we may therefore choose 
any number, jj {1 < j < N), such that 



(7) 



without any physical consequence. This ambiguity reflects the physical fact that a rigid stick 
responds identically regardless of how longifMdinal forces are distributed along the object. 
Inserting (6) and (7) into (1), we can thereby obtain the following system of equations 
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with 3(iV— 2)+l Lagrange multipliers, determined by the geometry of (2) and (3) (notice that 
the vector Afci = A^^ + A|j^ represents 3 scalar multipliers, one for the longitudinal direction 
A|i and two for orthogonal directions Xki)- These equations can now be introduced into 
a variety of numerical integrators, determining the Lagrange multipliers for enforcing the 
constraints, and exploiting the free Gauge parameters, 7j, for computational optimization. 



III. USING THE VELOCITY EXPLICIT VERLET METHOD 

In the following we illustrate the approach through the commonly used [2,4] velocity 
Explicit Verlet algorithm (shown here including linear damping) which approximates the 
solution to (1) using a non-zero time step of dt connecting time tn — ndt to tn+i — {n+l)dt, 

^r-Ri+{^-'^)dtvr+^^^ (9) 

= - 2_yn , ^ ^ ( j:n , pi+l\ QQN 

i -|- 2 2 ^ 

where we use the notation, X'^ — Xj{tn) — Xj{n • dt), to describe the integer n time step, 
and where Vj — Rj. 
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A. Determining the Lagrange Multipliers \ji 

Inserting (8) into (9) yields the constrained discrete-time equations 
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with = ^ and, following the notation of [2,4], R] is the updated (vector) position 
of particle j, had there been no constraints. The longitudinal and orthogonal components 
of Xji are here relative to ^i^v- The Lagrange multipliers are now determined such that the 
constraints are satisfied at time tn+i] i.e., under the conditions 
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Inserting (11) into (12) provides us with one equation for the determination of Aijv, 
p 
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whereas inserting (11) into (13) provides us with 3(A'" — 2) equations for the determination 
of 
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where 
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Solving equations (14)-(16) will provide the Lagrange multipliers. However, since (14) is a 
nonlinear equation, and since (14)- (16) are coupled, the direct solution to the equations will 
require an iterative approach to the solution of the 3(A^ — 2) -|- 1 equations. A significant 
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computational consequence of separating the orthogonal and longitudinal components of the 
Lagrange multiplier \ji is that the "free" parameters 7j (for 1 < j < N) introduced in (7) 
can now be chosen to optimize efficiency. Specifically, choosing jj = Ai/ An = rriN/mi leads 
to the following efficient, non-iterative, exact, and self- correcting procedure for evaluating 
the necessary multipliers: 

Step 1: Solve the two independent sets of — 2 linear equations (15) to obtain Xj^. 

Step 2: Solve (14). With A^^ evaluated, then p is known. Given = Ai/An, this equation 
becomes a simple second-order polynomial in Aijv where the root closer to zero is the 
relevant one. 

Step 3: With the completion of steps 1 and 2, A^^ can be found from (16) (with 7j = 
Ai/An) as the solution to the set of A'" — 2 linear equations. 

The facts that the method outlined above is non-iterative and exact are self-evident. Given 
that the involved matrices are rather small ((A^ — 2) x (A^ ~ 2)), the method is also com- 
putationally efficient. Finally, it can be seen that regardless of deviations from the desired 
constrained geometry at time n ■ dt, the geometric conditions for evaluating the Lagrange 
multiphers are enforced at time {n-\-l) ■ dt (see (12) and (13)). Thus, the method is self- 
correcting when making the important distinction between, e.g., a desired length. Iji and an 
actual length \R^i\, which may not be equal to Iji due to computational precision errors. 

The matrices corresponding to the N — 2 equations in (15) and (16) arc generally very 
well conditioned and are not subject to numerical instabifities. Specifically, for a system 
of A^ equal masses the matrix described in (16) has condition number [10] kl — (for 
A^ > 3). The matrix (15) describing the orthogonal components of the Lagrange multipliers 
is less transparent. However, for a system of evenly spaced equal masses, we show in Figure 
2 how the condition number /c^ depends on A^. It is clear from these data that the matrices 
are easily handled numerically for any reasonable number of constrained particles. It is 
important to realize that the matrix (15) may become less well conditioned if An/Ijn is 
extremely large for any < j < N. However, (11) and (15) have been expressed in terms 
of Xji, through the relationships (6) and (7). Thus, we could equally well have chosen the 
sought after Lagrange multipliers to be A^at, in which case equations equivalent to (15)-(17) 
can be produced to avoid a possible ill conditioned matrix for large An/Ijn- It is also 
important to realize that these matrices are not changing during the course of a simulation, 
and that, e.g., the same LU factorization can be effectively used for all time steps. The 
solution to (14) is limited mainly by p, as it can be seen that the polynomial may have 
only complex solutions when p is near orthogonal to Rin- For reasonable time steps, when 
particle positions only change by a fraction of particle size, p is always near parallel to Rin^ 
and (14) does therefore not pose any unpredictable numerical problems. 

The completion of the time step in (11) can then be accomplished for all participating 
particles in an ensemble. With these new coordinates at time ^n+i, the resulting forces f^'^^ 
can be evaluated. 
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B. Updating the Velocities 



As pointed out in [4], straightforward application of (10) will not in general satisfy the 
necessary constraints on the velocity complements to (12) and (13): 

(18) 
(19) 
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1< J < 



IN 



The necessary corrections can be introduced by the following modification, 
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where Bj — dt/{2mj + ajdt). As for the Lagrange multipliers Xji the introduction of cr^j 
should conserve angular momentum, i.e.. 
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where aji = aji + a^^ with (rfi-Rijq = and Ujj^\\RiN ■ Analogous to the physically arbitrary 
coefficient 7^, here we have coefficients 6j that can be chosen to be any number. With this 
observation, (20) becomes the velocity complement to (11): 
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The (vector) velocities without the corrective multipliers a are denoted VJ^~^^- The velocity 
equations are now determined by inserting (23) into (18) and (19), yielding the following set 
of 3(A" — 2) + 1 hnear equations. 
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where 



pn+1 



^ ^V''+^ - ^ ' ^^^--^v . (27) 



The solution to the set of hnear equations (24)-(26) can be significantly assisted by choos- 
ing 5k = Bi/Bn- This choice will allow (25) to become independent of the other equations, 
and all corrective multipliers a can be obtained by solving three sets of iV — 2 (plus one) 
linear equations in the hsted order. With this solution, the velocities VJ^'^^ can be completed 
through (23), and the time step is then accomplished. Before initiating another time step, it 
is important to include the corrective multipliers into the evaluated force, i.e., f^~^^ — ^j*"*"^, 
before evaluating i?""*"^. 

Notice that for small dt, the matrices given by (24) and (26) are nearly identical to the 
ones described by (15) and (16), respectively. Thus, we can assume that (24) and (26) are 
subject to very similar conditioning as discussed in the previous section, when solving the 
systems numerically. 



IV. CONCLUSION 



We have presented an efficient algorithm for conducting dynamical simulations of par- 
ticles constrained in a linear geometry. Following other presentations on geometrical con- 
straints in particle simulations, we introduce Lagrange multipliers to the equations of motion 
for each particle, and determine the multipliers by the desired constraints as well as basic 
conditions of conservation of momentum and angular momentum. The resulting equations 
for determining the Lagrange multipliers are, in general, a set of coupled equations with some 
nonlinear component, necessitating an iterative approach to obtain an approximate solution. 
By taking advantage of the physical ambiguity left in the distribution of the constraining 
forces in the longitudinal direction of the ensemble, we are able to determine a particular 
distribution that decouples the nonlincarity of the determining equations into a second-order 
polynomial of a single variable, namely, one of the Lagrange multipliers in the longitudinal 
direction. The rest of the equations are linear and can be solved as three independent sets 
of — 2 equations with N — 2 unknowns, where N is the number of mutually constrained 
particles. Thus, the algorithm is exact (within the accuracy of floating-point arithmetic), ef- 
ficient, and non- iterative. We have illustrated the approach through an implementation with 
the widely used velocity explicit Verlet algorithm for temporal integration of second-order 
differential equations. The presented formulation of the approach has the added benefit of 
being self-correcting, i.e., any spurious deviations in the constrained geometry arising from, 
e.g., computational precision error, is self-corrected, and the algorithm does not accumulate 
or propagate error in the constraints. It is straightforward to apply the method to other 
numerical methods for ordinary differential equations. For example, a recent comprehensive 
study on self-assembled monolayers [11] applied the above-outlined technique to model rigid 
and linear alkanethiol molecules as overdamped stochastically driven objects. 
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FIGURES 




FIG. 1. Sketch of the system under consideration. N particles are constrained to the hne, 
defined by Rn — Ri, through the Lagrange multiphers Xji = ^jl + ^ji, where 1 < j < N and 
i = l,N. 
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FIG. 2. Condition numbers, k2 and fe^, for the matrices described by equation (15) (o) and 
equation (16) (•) for 7^ = AyjA^ as a function of the matrix size {N — 2) x {N — 2), if the N 
particles have equal mass and are evenly spaced along Rin- The condition numbers are k2 
and k^p^{N- 2y-^. 
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